from numpy import floor, pi, sin, cos, degrees, radians, arcsin, arccos


class Solar_angle:
    def __init__(self, time_zone, year, month, day, hours, minutes, seconds, long, lat, Declination):
        self.time_zone = time_zone
        self.year = year
        self.month = month
        self.day = day
        self.lat = lat
        self.long = long
        self.hours = hours
        self.minutes = minutes
        self.seconds = seconds
        self.Declination = Declination

    def Calcu(self):
        day_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
        if self.year % 400 == 0 or (self.year % 4 == 0 and self.year % 100 != 0):
            day_month[1] = 29
            sum_day = sum(day_month[:self.month-1]) + self.day
        else:
            sum_day = sum(day_month[:(self.month - 1)]) + self.day

        beijing_hms = self.hours + self.minutes / 60 + self.seconds / 3600
        if self.long > 0:
            if self.time_zone == -13:
                long_cha = self.long - (floor((self.long * 10 - 75) / 150) + 1) * 15.0
            else:
                long_cha = self.long - self.time_zone * 15.0
        else:
            if self.time_zone == -13:
                long_cha = (floor((self.long * 10 - 75) / 150) + 1) * 15.0 - self.long
            else:
                long_cha = self.time_zone * 15.0 - self.long

        Ping_suntime = beijing_hms + long_cha / 15

        N_0 = 79.6764 + 0.2422 * (self.year - 1985) - int(0.25 * (self.year - 1985))
        t = sum_day - N_0
        theta = 2 * pi * t / 365.2422
        date_t = 0.0028 - 1.9857 * sin(theta) + \
                          9.9059 * sin(2 * theta) - \
                          7.0924 * cos(theta) - \
                          0.6882 * cos(2 * theta)

        truth_suntime = Ping_suntime + date_t / 60

        time_angle = (truth_suntime - 12) * 15

        if self.Declination is None:
            N = sum_day
            N0 = 79.6764 + 0.2422 * (self.year - 1985) - int((self.year - 1985) / 4)
            theta = 2 * 3.1415926 * (N - N0) / 365.2422
            Declination_angle = 0.3723 + 23.2567 * sin(theta) + 0.1149 * sin(2 * theta) - 0.1712 * sin(3 * theta) - \
                            0.758 * cos(theta) + 0.3656 * cos(2 * theta) + 0.0201 * cos(3 * theta)
        else:
            Declination_angle = self.Declination

        sin_alpha = sin(self.lat * pi/ 180) * sin(Declination_angle * pi/ 180) + \
                         cos(self.lat * pi/ 180) * cos(Declination_angle * pi/ 180) * \
                         cos(time_angle * pi/ 180)
        son_alpha = degrees(arcsin(sin_alpha))

        B1 = (sin(radians(Declination_angle)) * cos(radians(self.lat)) -
                   cos(radians(time_angle)) * cos(radians(Declination_angle)) * sin(radians(self.lat))) / cos(radians(son_alpha))
        B2 = (sin(radians(Declination_angle)) -
                   sin_alpha * sin(radians(self.lat))) / (cos(radians(son_alpha)) * cos(radians(self.lat)))

        if time_angle < 0:
            Solar_Azimuth = degrees(arccos(B2))
        else:
            Solar_Azimuth = 360 - degrees(arccos(B2))

        return time_angle, Declination_angle, son_alpha, Solar_Azimuth


if __name__ == '__main__':
    Solar_angle = Solar_angle(8, 2022, 7, 28, 10, 24, 10, 116.3833333, 39.9, Declination=18.932188)
    time_angle, Declination_angle, son_alpha, Solar_Azimuth = Solar_angle.Calcu()
